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This  annual  report  summarizes  progress  during  the  past  year. 


AERO-OPTICS  ANALYSIS 


Contract  DASG60-89-C-0145 
INTRODUCTION 

The  U.s.  Army  Strategic  Defense  Command  (USASDC)  has  several 
ongoing  and  planned  programs  that  utilize  optical  sensors  aboard 
missiles  traveling  at  hypersonic  velocities  in  the  atmosphere. 
Central  to  the  missile  homing  problem  are  aero-optical  effects 
upon  a  missile-borne  sensor/seeker  which  looks  through  both  an 
electromagnetic  window  and  the  flow  field  about  the  vehicle. 
Aspects  of  the  problem  include  modeling  and  simulation  of  the 
flow  field,  modeling  and  simulation  of  the  effects  of  the  flow 
field  on  incident  radiation  from  a  target,  and  finally, 
predicting  the  resultant  image  imperfections  and  error  in 
apparent  object  position  as  perceived  by  the  sensor. 

The  University  of  Alabama  in  Huntsville  is  engaged  in  a  program 
of  aero-optics  research  and  development  in  which  the  broad 
objectives  are  to:  develop  and  maintain  an  aero-optics  data 
base?  investigate  key  issues  in  the  areas  of  aero-optic  effects; 
and  increase  synergism  among  ongoing  efforts  in  aero-optics  both 
nationally  and  internationally  if  applicable. 

The  aero-optics  data  base  is  a  computerized  data  management 
system  to  house  aerodynamic  data  related  to  optics,  optics  data 
related  to  aerodynamics,  and  aero-optics  interrelated  data  that 
is  applicable  to  USASDC  projects  and  missions.  Data  shall  be 
from  both  test  and  analysis.  UAH  will  collect  and  update 
pertinent,  existing  data  and  will  update  the  data  base  with 
appropriate  new  data  as  these  data  are  obtained/generated. 
Further,  UAH  will  develop  a  data  acquisition  format  for 
laboratory  tests,  wind  tunnel  tests,  sled  tests  and  flight  tests 
to  meet  the  needs  of  potential  users,  e.g. ,  sensor  designers 
and/or  sensor  requirements  generators.  Moreover,  UAH  will  also 
develop  criteria  for  the  admission  of  data  to  the  data  base. 

Support  for  data  base  development  will  come  from  three  other 
phases  of  the  program,  namely,  analytical  investigations, 
assessment  of  measurement  techniques,  and  experiments  on 
hypersonic  shear  flow  simulations.  As  part  of  the  analytical 
effort,  UAH  is  surveying,  reviewing,  collecting,  updating  and 
using  available  aero-optical  mathematical  models  and  computer 
simulations.  In  assessing  measurement  techniques,  UAH  is 
concentrating  its  effort  on  non-intrusive  schemes  such  as 
holography,  optical  tomography,  laser  induced  fluorescence  and 
improved  standard  flow  visualization  such  as  Schlieren  and 
shadowgraph.  Finally  in  the  experimental  support,  UAH  is  setting 
up  a  laboratory  experiment  which  will  utilize  a  flow  channel 
having  low  supersonic  capability.  This  will  be  used  for  studies 
and  assessment  of  non-intrusive  measurement  techniques. 


The  following  sections  of  this  report  have  been  prepared 
by  individual  investigators  and  convey  progress  on  various  facets 
of  the  overall  program.  Each  section  is  self-contained  and  can 
be  read  as  an  entity  in  itself.  Specifically,  these  are 
presented  as  follows: 


*  Aero-Optics  Analytical  Studies  -  Dr.  c.  P.  Chen 

*  Hon-lntrusive  Measurement  Techniques  -  Dr.  J.  E.  smith,  Jr. 

*  UAH  Aero-optics  Data  Base  -  Dr.  J.  Ziebarth 

Dr.  S.  Graves 

*  optical  Measurement  Techniques  -  Dr.  M.  Abushagur 

*  Image  Restoration  by  Least  Squares  (Wiener)  Filter  - 
Mr.  K.  S.  Rim 

*  Mutual  Coherence  Propagation  in  Turbulent  Media  - 
Mr.  A.  Monteiro 

Dr.  J.  Jarem 


Aero-Optics  Analytical  Studies 
by 

c.  P.  Chen 


Objective 

The  objective  of  the  analytical  studies  is  to  survey  and  evaluate 
current  computational  fluid  dynamics  codes  for  predicting 
turbulent  shear  layers  above  optical  windows.  Current  state-of- 
the-art  physical  modeling,  in  terms  of  turbulence  and  real  gas 
effects,  is  also  reviewed  and  assessed. 

Technical  Problems 

The  turbulent  mixing  shear  layers  over  optical  windows  (created 
by  a  slot  wall  cooling  jet)  of  hypersonic  vehicles  cause  severe 
limitations  to  target  detection  and  recognition.  Further,  if  the 
wall  jet  nozzle  is  not  properly  expanded,  a  series  of  oblique 
shocks  and  expansion  wave  can  also  form  above  the  window.  Thus, 
in  a  coupled  aerodynamic-optical  analysis,  the  fluid  dynamics 
analysis  must  predict  not  just  pressure  and  velocities,  but  also 
steady  and  unsteady  density  fields.  In  aero-optics  applications, 
variable  density  effects  are  primarily  caused  by  the  mixture  of 
gases  with  different  densities,  strong  temperature  gradients 
within  the  fluid,  strong  distortion  of  shock-boundary  layer 
interactions  and  compressibility  effects  in  high  speed  flows. 
The  Mach  number  and  Reynolds  number  of  a  typical  operating  flight 
envelope  precluded  the  use  of  most  traditional  aero-dynamic 
design  methodologies.  In  such  approaches,  simplified  analyses 
are  combined  with  experimentally  generated  engineering 
correlations.  To  supplement  or  replace  these  traditional  design 
methods,  computational  fluid  dynamics  (CFD)  techniques  are  being 
utilized  as  a  design  tool  over  portions  of  the  flight  envelope 
where  experimental  data  cannot  be  obtained  economically  or 
feasibly.  Despite  the  current  maturity  of  state-of-the-art  CFD 
codes,  much  work  remains  to  be  done  before  they  can  be  used  with 
confidence  in  a  design  or  diagnostic  environment. 

General  Methodology 

The  first  step  in  applying  CFD  methods  to  Aero-optics  flow 
applications  was  to  conduct  a  survey  to  identify  CFD  codes  in  the 
public  domain,  i.e.,  available  from  government,  industries  and 
university  sources.  Due  to  strong  shear  and  turbulence  effects 
in  the  flow  fields  related  to  aero-optics,  Euler  codes  such  as 
CM3D  and  Eagle  were  not  included.  The  capabilities  of  the 
applicable  viscous  codes  identified  in  this  survey  are  summarized 
in  Table  1.  In  selecting  these  codes,  special  attention  was  paid 
to  the  code's  capabilities  in  handing  turbulence  and  real  gas, 
chemical  reacting  flows.  Both  Parabolized  Navier-Stokes  (PNS) 
codes  and  full  Navier-Stokes  (NS)  codes  were  surveyed.  PNS 
techniques  can  be  applied  to  both  the  invicid  outer  flow  as  well 
as  the  viscous  inner  flow.  Although  such  methods  are  still  not 


applicable  to  flow  with  strong  interaction  effects  such  as 
separation,  they  are  ideal  analytical  tools  for  their  efficiency 
and  modest  computer  requirements.  For  shear  layer  driven,  strong 
shock-boundary  layer  interaction  flows,  full  Navier-Stokes 
analyses  are  required. 

A  correct  description  of  the  turbulence  structure,  which  includes 
the  variable  density  efforts,  is  the  most  important  feature  of 
the  aero-optics  flow  field  analysis.  Current  state-of-the-art  in 
turbulence  modeling  cannot  generate  accurate  predictions  under 
general  conditions.  What  can  be  expected  is  that  satisfactory 
predictions  can  be  made  for  categories  of  flows  for  parameters  of 
aero-optics  concern,  even  if  that  requires  zonal  or  ad  hoc 
modeling  modifications,  so  long  as  the  model  predictions  are 
verified  as  acceptable.  Thus,  by  surveying  these  CFD  codes, 
advances  in  turbulent  fluctuation  predictions,  compressibility, 
and  turbulence/chemistry  coupling  can  be  recommended. 

Technical  Results 


For  the  CFD  codes  identified,  a  systematic  summary  of 
characteristics  in  terms  of  numerics,  physical  models  (turbulence 
models  and  chemistry  models)  and  validation  has  been  generated. 
For  many  CFD  codes,  a  sequential  development  in  terms  of  numerics 
upgrades  and  thermo-physical  model  upgraded  can  be  identified. 
In  such  cases,  hierarchial  characterization  approaches  are  also 
used.  This  summary  is  directly  communicated  as  input  to  the  Data 
Base. 

To  gain  some  confidence  in  using  CFD  codes,  validations  or 
calibrations  of  the  codes  against  benchmark  experimental  data  are 
necessary  to  examine  the  accuracy  of  the  physical  and  numerical 
models,  as  well  as  to  highlight  the  procedures  to  obtain  reliable 
results.  Due  to  the  lack  of  a  sufficiently  detailed  data  base 
describing  various  flow  processes  for  a  hypersonic  flight 
vehicle,  several  "unit  problems"  that  delineate  dominant  physical 
phenomena  applicable  for  specific  components  and  portions  of 
aero-optics  flows  are  identified.  These  include:  supersonic 
mixing  layer,  cavity  flows,  wall  jets  and  shock/boundary  layer 
interaction.  Comparison  calculations  using  different  codes 
versus  available  experimental  data  or  among  various  CFD  results 
are  then  made.  Code  to  code  comparisons  help  quantify  numerical 
errors  between  algorithms  when  identical  physical  models  are 
solved  with  different  methods.  The  series  of  validation 
calculations  on  chosen  unit  problems  in  terms  of  accuracy  and 
efficiency  are  also  summarized  as  inputs  to  the  Data  Base. 

Important  Findings  and  Conclusion 

In  examining  features  of  these  CFD  codes,  it  can  be  found  that 
each  code  assumes  some  degree  of  physical  or  numerical 
approximations.  Numerically,  most  codes  are  equipped  with  shock 
capturing  capabilities  by  using  flux-splitting,  TVD  or  some 
dissipation  schemes.  Most  codes  also  use  the  so-called  density- 
based  method  in  which  the  governing  equations  are  solved  in  a 
coupled  fashion.  However,  a  majority  of  these  codes  assume  that 


the  flow  obeys  a  perfect  gas  low  and  also  assume  calorically 
perfect  thermodynamics.  Only  a  limited  number  consider  the 
effects  of  air  chemistry  on  the  gas  properties.  However,  the 
typical  flow  environment  is  such  that  ideal  gas  assumptions 
cannot  be  applied  without  careful  consideration  of  the  flight 
envelope  and  the  specific  geometry. 

The  major  shortcomings  of  the  most  of  these  codes  are  the 
turbulence  models  used  in  CFD  predictions.  Turbulence  in  aero- 
optics  flows  calls  for  predictions  requirements  for  both  the  mean 
flow  and  fluctuating  velocity  and  density  field.  All  the  codes 
use  one-point  closure  methods  which  include  algebraic  models 
and  one,  as  well  as  two  equation  models.  Turbulence/chemistry 
interactions  are  not  accounted  for  in  all  codes.  Unsteady 
effects  are  also  not  modeled  properly. 

To  conclude,  application  of  CFD  methods  to  analyze  aero-optics 
flows  is  a  powerful  tool  for  supplementing  and  extending 
traditional  design  methods.  However,  the  proper  use  of  CFD 
requires  a  systematic  application  methodology  that  uses  existing 
codes  where  possible  and  includes  improved  numerics  and  physical 
models  where  dictated  by  the  geometry  and  flow  regimes  under 
consideration.  These  codes  should  then  be  exercised  using  a 
series  of  validation  calculations  on  appropriately  chosen  unit 
problems  to  determine  the  applicability,  accuracy  and  efficiency 
of  the  code. 


< 

r  ^ 

•A 

cn 

Cu 

£ 

*  % 
s-/ 

'w 
»— • 

cn 

2 

cn 

Ul  Ci. 

cz 

as 

cH 

2 

Q 

<$ 

51  ** 

< 

< 

< 

< 

o 

O 

r. 

O 

<1 

21  - 

ri 

C~ 

-r 

lO 

»A 

**0 

co 

ci 

o 


turbulence  model  for  hypersonic  How  calculations. 


Objective 


Non-intrusive  Measurement  Techniques 

by 

J.  E.  Smith,  Jr. 


The  objective  of  this  phase  of  the  program  is  to  explore  and/or 
develop  non-intrusive  measurement  techniques  under  high  speed 
flow  conditions  to  characterize  wall  jet  and  mixing  layer 
phenomena  which  influence  the  quality  of  aero-optic  measurements. 


Approach 

The  basic  approach  consists  of  two  principal  elements  as  follows: 

(  i)  Develop  a  computer  controlled  experimental  apparatus  to 
perform  short  duration  wall  jet  experiments. 

(ii)  Examine  the  flowing  stream(s)  with  various  non-intrusive 
measurement  techniques  to  evaluate  local  speciation  and 
velocity.  Initially,  we  will  investigate  the  use  of 
chemical  tracers  and  Raman  Spectroscopy  conducted  over  fiber 
optics  to  determine  local  speciation  within  the  wall  jet. 


Status 

A  schematic  of  the  flow  system  is  given  in  the  accompanying 

figure.  Relative  to  that  sketch,  the  brief  status  is  as  follows: 

•  Pressure  supply  system,  vacuum  system  and  scrubber  are  90% 
complete. 

*  The  test  section  is  under  design. 

•  Spectroscopic  unit  is  75%  complete. 

*  The  LDV  unit  will  require  software  and  hardware  upgrades  to 
make  it  fully  compatible  with  the  existing  data  acquisition 
and  computer  equipment. 


Future  Plans 

These  plans  include  the  following  elements: 

"  Complete  construction  of  flow  system  and  sub-systems. 

•  Complete  CAD  design  of  the  test  section  and  manufacture  the 
design  using  a  computer  driven  milling  machine  available  at 
the  university. 

•  Conduct  Raman  Scattering  experiments  to  measure  the  local 
concentration  of  a  chemical  tracer  added  to  the  wall  jet. 

•  Upgrade  the  LDV  unit  to  evaluate  hybrid  system  to  measure 
local  velocities  with  LDV  and  concentrations  with  Raman 
Spectroscopy. 
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Objective 


UAH  Aero-Optics  Data  Base 
by 

J.  Ziebarth  and  S.  Graves 


The  basic  objective  is  to  develop  a  computerized  data  management 
system  that  will  house  aerodynamic  data  related  to  optics  and 
optics  data  related  to  aerodynamics. 

Approach 

In  order  to  accomplish  this  task,  the  fundamental  approach  is 
given  as  follows: 

’  Integrate  a  DBMS  and  a  menu  driven  interface  system,  to 
provide  a  user  friendly  AODB  which  is  flexible  and 
extendible. 

’  Include  analytical/experimental  data,  visualization/analysis 

tools,  bibliographical  information. 

’  Maintain  compatibility  with  USASDC-KDEC. 

Status 


The  prototype  data  base  is  complete  and  contains  representative 
data  to  demonstrate  functionality.  Moreover,  a  Configuration 
Board  has  been  established  to  direct  the  selection  criteria  and 
management  of  data. 

Future  Plans 

These  plans  consist  of  four  key  elements  as  follows: 

(  i)  Expand  capability  as  required  by  users. 

(  ii)  Acquire  HW/SW  for  a  fully  operational  AODB. 

(iii)  Provide  for  acquisition  of  existing  and  future  data. 

(  iv)  Establish  connectivity  between  UAH-AODB  and  USASDC-KDEC. 


Optical  Measurement  Techniques 
for 

Turbulent  Flow  Fields 
by 

M.  Abushagur 


The  turbulent  flow  in  front  of  the  window  of  an  interceptor 
consists  of  a  medium  with  variable  density  and  in  turn  with 
variable  index  of  refraction.  A  light  beam  wavefront  traveling 
through  this  medium  suffers  a  phase  change  that  deforms  the 
incoming  wavefront.  The  deformation  of  the  wavefront  affects  the 
quality  of  the  image  observed  through  this  medium.  The  poor 
quality  of  the  image  of  the  target  effects  the  tracking 
capability  of  the  interceptor.  Overcoming  such  a  problem 
requires  a  knowledge  of  the  flow  field  density  distribution  in 
front  of  the  window.  This  density  distribution  can  be  used  to 
compute  the  transfer  function  of  the  turbulent  medium.  The 
transfer  function  then  can  be  used  in  designing  the  proper 
imaging  enhancement  system. 

Overcoming  such  a  problem  requires  a  knowledge  of  the  flow  field 
density  distribution  in  front  of  the  window.  This  density 
distribution  can  be  used  to  compute  the  transfer  function  of  the 
turbulent  medium.  The  transfer  function  then  can  be  used  in 
designing  the  proper  imaging  enhancement  system. 

Measuring  the  density  distribution  without  perturbing  the  flow  is 
a  very  critical  issue  in  this  process.  Optical,  non-intrusive 
methods  are  adequate  to  such  measurements.  These  techniques  are 
based  on  the  fact  that  the  wavefront  of  the  light  beam,  as  it 
passes  through  the  turbulent  medium,  is  modulated  by  the  index  of 
refraction  distribution.  When  extracted,  this  information  can  be 
used  in  mapping  the  density  distribution  of  the  flow  field. 

There  are  a  number  of  optical  methods  which  can  be  used  for 
mapping  the  flow  field.  Among  these  methods  are: 
interferometry,  holographic  interferometry,  speckle  photography, 
shadowgraphs,  Schlieren  technique,  and  Moire  deflectometry. 
These  methods  are  not  equivalent,  each  has  its  own  advantages  in 
measuring  a  certain  quantity  of  the  flow  field;  they  also  differ 
very  significantly  in  their  complexity.  The  method  we  are 
interested  in  using  in  this  preliminary  stage  of  our  study  is 
Moire  deflectometry,  which  is  very  powerful  in  mapping  phase 
objects  such  as  the  case  for  flow  fields. 

(a)  TASK  OBJECTIVES 

The  objective  of  this  work  will  be  to  develop  optical  techniques 
to  reconstruct  the  flow  field  distribution  in  front  of  the 
interceptor  window.  To  achieve  this  objective,  a  number  of 
simpler  tasks  are  set  to  develop  the  measurement  techniques. 


These  subtasks  are  as  follows: 


1.  Develop  non-intrusive  optical  methods  for  measuring  simple 
flow  fields  which  are  known  and  for  which  theoretical  models 
have  been  developed  such  as:  the  temperature  distribution 
in  a  flame;  temperature  field  around  a  heated  wire;  two 
heated  wires ,  ...  etc . 

2.  Develop  the  software  capable  for  interpreting  the  optical 
data  collected  from  the  experiments  in  terms  of  flow  field 
distribution  plots. 

3.  Test  the  algorithms  developed  in  subtask  (2)  for  the 
measurement  of  the  flow  fields  of  the  models  mentioned  in 
subtask  (1) . 

4.  Build  a  system  (both  hardware  and  software)  capable  of 
measuring  flow  field  distribution  automatically,  i.e. 
collecting  data,  analyzing  it,  and  plotting  the  flow  field 
distribution. 

(b)  TECHNICAL  PROBLEMS 

The  flow  field  distribution  measurement  in  front  of  the 
interceptor's  window  is  not  feasible  practically.  Wind  tunnel 
measurements  are  the  closest  to  the  real  life  situation.  At  this 
stage  of  this  research,  we  are  considering  more  fundamental 
problems,  namely,  the  measurement  of  a  symmetric  flow  field 
distribution,  an  asymmetric  flow  and  finally,  a  mixing-layer  flow 
field. 

(C)  GENERAL  METHODOLOGY 

The  Moire  effect  denotes  a  fringe  pattern  formed  by  the 
superposition  of  two  grid  structures  of  similar  period  with  some 
inclination  angle  as  shown  in  Fig.  1.  This  method  is  used  in 
many  applications.  Yet  another  more  powerful  technique  for  the 
class  of  problems  we  are  considering  is  Moire  deflectometry.  In 
the  classical  Moire  fringes,  the  two  gratings  are  superimposed  on 
one  another.  On  the  other  hand,  in  Moire  deflectometry,  the  two 
grating  are  placed  apart.  In  this  technique,  the  deformation  of 
the  Moire  pattern  is  caused  by  the  ray  deflections  due  to  their 
passage  through  the  phase  object  (e.g.  flows  fields) .  In  this 
technique,  the  object  to  be  tested  is  placed  in  the  path  of  a 
collimated  beam  followed  by  a  pair  of  transmission  gratings 
(Ronchi  rulings)  placed  at  a  distance  from  each  other  as  shown  in 
Fig.  2.  The  Moire  deflectogram  (the  resulting  fringe  pattern)  is 
a  map  of  ray  deflection  corresponding  to  the  optical  properties 
(index  of  refraction  distribution)  of  the  test  object. 

The  Moire  fringe  pattern  is  produced  by  crossing  two  Ronchi 
rulings  as  shown  in  Fig.  1.  The  fringe  spacing  depends  on  the 
angle  8. 


p 


(1) 


sine 

where  p  and  p1  are  the  pitches  of  the  Ronchi  ruling  and  Moire 
fringe  spacing,  respectively.  These  fringes  can  be  seen  by 
illuminating  the  Ronchi  ruling  with  a  coherent  plane  wave.  Due 
to  the  Talbot  effect,  these  fringes  can  be  seen  at  periodic 
distances . 

If  the  plane  wave  is  obstructed  with  an  object  that  has  a 
variable  optical  pathlength  (optical  pathlength  is  the  product  of 
the  object  thickness  and  its  index  of  refraction)  the  fringes  no 
longer  are  straight.  The  shape  of  the  produced  fringes  depends 
on  the  index  of  refraction  distribution.  As  shown  in  Fig.  3, 
when  a  light  ray  passes  through  the  phase  object,  because  of  the 
index  of  refraction,  it  gets  refracted  (bent)  with  an  angle  d. 
This  angle  is  related  to  the  index  of  refraction  by  the  following 
equation: 


d(y)  = 


|  -2  n  (x,y ,  z) 


where  NQ  is  the  index  of  refraction  of  the  outside  medium,  and 
^  /  ”3  x1  is  the  partial  derivative  with  respect  to  x1. 


N(Y,  X)  =  n(x,y, z)  -NQ, 
Then  N(Y,  X)  can  be  shown  to  be 
N(Y,X)  =  — L  jv*  „Q  f 


dtY1) 

rsin (x-e) -y1 


dy1. 


Hence,  the  index  of  refraction  distribution  N(r,x)  can  be 
computed  from  Eg.  (4)  using  d(yi)  data.  The  deflection  angle 
distribution  d(y  ,  z)  can  be  determined  from  the  intensity 
distribution,  I(y,z),  of  Moire  deflectogram  from  the  following 
relation. 

d (y1 ,  z )  =  mirl  (y^z)  +  2k7r?  (5) 

where  m  =  II  and  k  *  0,  II,  12,  ... 

The  values  for  m  and  k  can  be  determined  by  a  set  rules.  After 
determining  the  deflection  angle  distribution,  the  index  of 
refraction  distribution  can  be  calculated  from  the  inverse  of  the 
Abel's  transform,  Eq.  (4). 

(d)  TECHNICAL  RESULTS 


An  experiment  was  performed  at  UAH  ECE  labs  for  measuring  the 
temperature  distribution  across  a  flame.  The  setup  used  is  shown 
in  Fig.  4.  A  HeNe  laser  was  used  to  illuminate  a  candle  flame; 
the  light  is  then  passes  through  two  Ronchi  rulings  spaced  at 
their  first  periodic  distance.  A  photodiode  was  used  to  measure 


the  light  intensity  across  the  Moire  deflectogram.  The  intensity 
distribution  is  plotted  in  Fig.  5  as  a  function  of  the  distance 
across  the  deflectogram.  This  intensity  distribution  data  was 
then  used  to  determine  dfy^z)  in  Eq.  (5).  The  deflection  angle 
computed  is  plotted  in  Fig.  6.  The  data  from  the  deflection 
angle  is  substituted  in  Eq.  (4)  to  determine  the  inverse  Abel's 
transform  to  compute  the  index  of  refraction  distribution  through 
the  flame;  this  distribution  is  inversely  proportional  to  the 
temperature  distribution.  The  temperature  distribution  which 
results  from  these  calculations  is  shown  in  Fig.  7. 

The  resultant  temperature  distribution  across  the  flame  is  a  good 
representation  and  agrees  with  theoretical  and  direct 
measurements . 

(e)  IMPORTANT  FINDINGS 

The  Moire  deflectometry  method,  used  in  this  stage  of  research 
for  the  flow  field  distribution  measurement,  was  used  to 
determine  the  flow  field  (temperation  distribution)  of  a  flame. 
Results  were  very  encouraging.  This  method,  in  the  future,  will 
be  extended  to  different  objects  as  outlined  previously  in 
section  (a)  of  this  report. 
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PREFACE 


This  document  is  a  report  on  restoring  image  which  is  biurred  by  Aero-Optics  effects  via  Least- 
Squares  (Wiener)  Filter.  The  goal  of  this  research  is  to  develope  a  digital  filter  algorithm  for 
<\stimating  the  ideal  image  from  the  blurred  image  using  a  posteriori  information.  Technical 
work  was  accomplished  by  Mr.  Kveong-Seop  Kim.  a  graduate  student  in  the  ECE  Department, 
under  the  supervision  of  Dr.  Reza  Adhami. 
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1  Introduction 


Images  are  produced  to  record  or  display  useful  information,  but  the  process  of  image  formation 
and  recording  is  imperfect.  The  recorded  image  invariably  represents  a  degraded  version  of 
the  original  scene.  Blurring  can  be  caused  by  relative  motion  between  the  camera  and  the 
original  scene,  or  by  an  optical  system  that  is  out  of  focus.  When  aerial  photographs  are 
produced,  blurs  are  introduced  by  atmospheric  turbulence,  and  relative  motion  between  the 
••amera  and  the  object.  An  example  of  such  blurring  process  is  the  Aero-Optic  Mixing/ Shear 
Layer  measurements  problem.  When  a  hypersonic  vehicle  is  launched,  the  window  in  this  vecjt^lg 
must  be  cooled.  Unfortunately,  the  mixing/shear  layer  formed  between  the  environment  gas 
and  the  cooiant  gas  causes  severe  optical  degradation.  Thus  we  need  to  restore  or  estimate 
the  original  scene  from  this  degraded  one  by  gas  turbulence.  The  field  of  image  restoration  is 
•■oncerned  with  the  reconstruction  or  estimation  of  an  uncorrupted  image  from  a  distorted  and 
noisy  one.  Thus,  restoration  techniques  are  oriented  toward  modeling  the  the  degradation  and 
applying  the  inverse  process  in  order  to  recover  the  original  image.  This  report  is  arranged 
into  several  sections.  Section  2  discusses  mathematical  models  for  images  and  blur  operators. 
Procedures  for  deblurring  require  complete  knowledge  of  the  blurring  function.  As  this  is  rarely 
available.  Section  3  discusses  parametric  blur  identification  for  the  estimating  of  the  blurring 
operator  from  the  blurred  image  itself.  Section  4  discusses  the  Wiener  image  restoration  filter. 
The  Wiener  filter  is  the  optimum  filter  for  minimizing  errors  between  the  given  degraded 
image  and  the  estimated  ideal  image.  The  data  obtained,  the  experiments  performed  including 
nonparametric  blur  identification  and  Wiener  filtered  image  restoration,  and  the  conclusions 
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drawn  from  this  research  are  listed  in  Section  5. 


2  Models  for  Blurred  Image  Formation 

2.1  Image  Formation 

For  notational  convenience,  consistency,  and  ease  of  understanding,  the  following  format  will 
he  established:  * 

•  g(x,y)  will  be  the  recorded  image; 

•  /(x,  ?/ )  will  be  the  ideal  image,  which  is  a  2-D  mapping  of  the  3-D  input  scene; 

•  h(x.y,s,t)  will  be  the  2-D  impulse  response  (point  spread  function): 

•  .?{.}  will  be  the  sensor  nonlinearity  which  wiil  be  modeled  as  a  point  operator: 

•  nix.y)  will  be  2-D  noisy  process: 

If  the  image  formation  process  is  linear,  the  recorded  image  can  be  modeled  as  the  output  of 
the  system  shown  in  Fig.  i,  which  is  given  mathematically  by 

f-roo  r+co 

g{x,y)  =  s{  /  h(x, y\s,  t)f(s.  t)dsdt)  +  n(x.y)  (1) 

J—  00  j  — OO 
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■f(x,y) 


n(x,y)  Noise 


9(x,y) 


Image  Blurring 


Sensor  Blurred  Image 

Figure  i:  Model  for  the  processes  of  image  formation  and  recording. 

The  noise  contribution  is  shown  in  Fig.l  as  an  additive  random  process  which  is  stalls- 
ticaliy  uncorrelated  with  the  image.  If  the  impulse  response  is  stationary  across  the  image  and 
object  fields,  it  becomes  a  function  of  only  the  argument  differences  z  -  s  and  y  -  t.  In  this 


case. 


(2) 

(3) 


9iz’Vl  =  5</_„  — <)/( a,t)dsdt}  +  n(i,y) 

g{x,y)  =  *{h(x,y)  *  f(x,y)}  +  n(x,y) 
where  (")  is  used  to  denote  2-D  convolution. 

in  a  discrete  impiementation.  continous  arguments  f.  g,  h  and  n  are  replaced  by  arrays  of 

-ampies  taken  on  N  <  .V  2-D  rectangular  lattices  of  equi-spaced  samples.  The  sampled 
are  related  bv 


arrays 


.v-i  ,v-i 

9{iJ)  £  Mi'j]kJ)f(k,l)}  +  nU',;') 


0  <  ij  <  .V  -  1 


4=0  (=0  "  "  '  -  -  j  »  (4) 

For  the  spatially  invariant  (stationary)  system,  the  convolution  integral  (2)  becomes  a 
convolution  sum 


.v-i  ,v-i 

s(i\j)  =  H  £  £  *(<  -  k.j  -  l)f(k.l))  +  n (i.i) 

4=0  (=0 

g(ij)  =  -s  WJ)*f(iJ)}  +  n(i,j) 


(5) 

(6) 
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where  the  asterisk  (*)  is  now  used  to  denote  a  discrete  convolution.  Often  the  sensor  nonlinearity 
is  conveniently  neglected  (or  linearized)  to  justify  the  use  of  a  linear  restoration  filter.  When 
this  nonlinearity  is  ignored.  (6)  reduces  to  the  linear  convolution  model 

g{hi)  =  h{i,j)  *  f(i,j )  +  n{i,  j)  (7) 

for  which  discrete  Fourier  transform  can  be  used  to  yield  the  frequency  domain  model 

G(m.n)  =  H(m.n)F(m.n i  -f  N(m.n)  (8) 

Here  H(m.n)  represents  samples  of  the  frequency  response  of  the  blurring  system  and  m  and 
n  are  the  discrete  frequency  variables. 

2.2  Image  Model 

Certain  linear  image  restoration  techniques  including  Wiener  filter  make  use  of  a  priori  sta¬ 
tistical  knowledge  of  the  original  (undistorted)  image.  This  takes  the  form  of  a  power  density 
spectrum  for  the  Wiener  filter.  A  large  class  of  real-world  images  can  be  modeled  as  the 
following  2-D  autoregressive  process  of  low  order: 

/(*.;)=  £  «(p, <?)/(» -/'J  -*/)  +  «(*, j)  .  V(i, j)  (9) 

Here  u(i,j)  can  be  viewed  as  the  error  in  approximating  f{i,j)  using  a  linear  combination  of 
neighboring  sample  values  contained  in  a  neighborhood  IFi.  Some  common  choices  for  W\  are 
illustrated  in  Fig.' 2.  A  comprehensive  survey  of  these  three  image  models  has  been  given  by 

[5]- 
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(a) 


(c) 


(b) 

Figure  2:  Image  model  support  corresponding  to  (a)  nonsymmetric  halfplane:(b)  sexnicausal;(c) 
noncausal  image  model. 

2.3  Blur  Models 


•  Motion  Blur:  Many  types  of  motion  blur  are  caused  by  relative  motion  between  the 
camera  and  the  object.  This  can  be  in  the  form  of  a  translation,  a  rotation,  a  sudden 
change  of  scale,  or  some  combination  of  these.  For  example,  when  the  object  translates 
at  a  constant  horizontal  velocity  V  during  the  exposure  interval  [0,T],  the  distortion  is 
one-dimensional  and  the  discrete  point-spread  function  makes  use  of  the  blurring  distance 
L,  which  is  the  number  of  additional  points  in  the  image  resulting  from  a  single  point  in 
the  original  scene. 

HhJik'l)  =  77~T\  •  0  <i-k<L  (10) 

\  T  I ) 


The  frequency  response  corresponding  to  this  blur  is  given  by 


1  ,  -,Lts 

H(m'ni = (in 


This  frequency  response  shows  zero  on  lines  parallel  to  the  n-axis  with  an  interline  spacing 
of  The  presense  of  these  parallel  zeros  in  the  frequency  domain  indicates  the 

direction  of  motion. 
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•  Out-of- Focus  Blur:  When  a  three-dimensional  scene  is  imaged  by  a  camera  on  some  parts 
of  the  scene  are  in  focus  while  other  parts  are  not.  The  degree  of  defocus  depends  upon 
the  effective  lense  diameter  and  the  distance  between  the  object  and  the  camera.  The 
corresponding  point  spread  function  is 

h(x,y)  =  ,  x-  +  i r  <r2  (12) 

■xr1 

where  r  is  the  radius  of  the  circle  of  confusion.  The  frequency  response  H{m.n )  is  given 
by 

H(m.n)  =  ^  —  .  p~  —  rn'  -f  n*  (13) 

rp 

where  J\  is  the  first-order  Bessel  function. 

•  Blur  due  to  atmospheric  turbulence 

h{x,y)  -  C  exp (14) 

3  Parametric  Blur  Identification 

If  the  type  of  degradation  is  known  such  as  linear  motion  blur  and  out-of-focus  blur,  the  blur 
can  be  identified  by  the  parametric  blur  description.  For  linear  motion  blur,  as  given  in  (10)  it 
is  only  necessary  to  estimate  the  direction  of  blur  and  the  blurring  distance.  Such  an  motion 
blur  is  often  sufficient  to  estimate  the  zero  patterns  in  the  frequency  response  from  which  one 
can  estimate  the  direction  of  motion  and  the  blurring  distance.  With  the  simplified  model  for 
an  out-of-focus  blur  in  Eq.  (12),  it  is  only  necessary  to  estimate  the  radius  of  the  circle  of 
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confusion.  In  this  case,  the  frequency  zero  patterns  can  be  used  to  estimate  the  radius  of  the 


circle  of  confusion. 

4  Wiener  filter  for  restoring  blurred  image 

4.1  Inverse  Filter  for  Image  Restoration 

Let  f.  g,  and  n  represents  XX  x  1  column  vectors  formed  by  stacking  the  rows  of  the  N  x  N 
and  n(i,j).  The  rirst  X  elements  of  f,  for  example,  are  the  elements  in  the  first 
row  of  the  next  X  elements  are  from  the  second  row,  and  so  forth  for  all  N  rows  of 

f{x,y).  Using  this  convention.  Eq.  (5)  can  be  expressed  in  the  following  vector-matrix  form: 

g  =  Hf  +  n  (15) 

where  H  is  of  dimension  A’2  x  A’2.  More  details  for  this  description  are  explained  in  [4].  From 
Eq.  (15).  the  noise  term  in  the  degradation  model  is  given  by 

n  =  g  -  Hf  (16) 

In  the  absense  of  any  knowledge  about  n.  a  meaningful  criterion  function  is  to  seek  an  f  such 
that  Hf  approximates  g  in  a  least-squares  sense  by  assuming  that  the  norm  of  the  noise  term 
is  as  small  as  possible.  In  other  words,  we  wish  to  find  an  f  such  that 

l|n|f-=||g-Hf  ||2  (17) 
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is  minimum,  where,  j|  .  |j  is  norm  operator.  From  Eq.  (17),  we  may  equivalently  view  this 
problem  as  one  of  minimizing  the  criterion  function 

7(f)  =  a  g  -  Hf  f  (18) 

with  respect  to  f. 

Minimization  of  Eq.  ( IS)  can  be  obtained  by  simply  differentiate  J  with  respect  to  f.  and 
'■et  equal  to  zero,  that  is. 

=  0  =  — 2H'(  g  -  Hf)  (19) 

di 

Here  H’  means  the  transpose  of  H.  Solving  Eq.  (IS)  for  f  yields 

f  =  H-1(H')~1H,g  =  H-1g  (20) 

4.2  Wiener  Filter  for  Image  Restoration 

Let  Q  be  a  linear  operator  on  f.  In  this  section  we  consider  the  least  squares  restoration  problem 
as  one  of  minimizing  functions  of  the  form  j|  Q  f  ||2,  subject  to  the  constraint  ||  g  -  H  f||2  = 
||  n  ||2.  This  approach  introduces  considerable  flexibility  in  the  restoration  process  because  it 
yields  different  solutions  for  different  choices  of  Q.  The  addition  of  an  equality  constraint  in 
the  minimization  problem  can  be  handled  by  using  the  method  of  Lagrange  multipliers.  The 
procedure  is  to  express  the  constraint  in  the  form  q(||  g  -  H  f||2  —  |(  n  ||2)  and  then  append  it 
to  the  function  II  Qf  II2-  In  other  words,  we  seek  an  f  which  minimizes  the  criterion  function. 

J(f)  =!|  Qf  f  +c(||  g  -  Hf  II2  -  [|  n  II2)  (21) 
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where  ct  is  a  constant  called  the  Lagrange  multiplier.  Differentiating  Eq.  (20)  with  respect  to 

A  * 

f  and  setting  the  result  equal  to  zero  yields  the  solution  f;  that  is 


f  =  (H/H  +  7Q'Q)"1H'g 


(22) 


where  7  =  -k 

9  Ol 

Let  Rf  and  Rn  be  the  correlation  matrices  of  f  and  n.  defined  respectively  by  the  equations 

Rf=£{ff'}  ,  Rn  =  £{rm'}  (23) 

where  E  {.}  denotes  the  expected  value  operation.  The  ijth  element  of  Rj  is  given  by  E  {/,-/,}, 
which  is  the  correlation  between  the  ith  and  the  jth  elements  of  f.  By  defining 

Q'Q  =  Rf'Rn  (24) 

and  substituting  this  expression  in  Eq.  (22)  we  obtain 

f  =  (H'H  +  7Rf-‘RnrlH'g  (25) 


In  frequency  domain,  Eq.  (25)  is  transformed  to 


F(u,u)  =  ( 


H‘{  u,u) 


•)G(it,u) 


(26) 


for  u.  v  =  0.  1,  2, .  Ar-1.  where  \H{u.v\2  =  Hm{u,v)H(u.v)  and  5/(u,u)  =  E{R/},  Sn(u,v) 

=  F{Rn}.  when  7=1,  the  term  inside  the  outer  brackets  in  Eq.  (26)  reduces  to  the  so-called 
Wiener  filter.  In  the  absense  of  noise  (5n(u,u)  =  0)  ,  the  Wiener  filter  reduces  to  the  ideal 
inverse  filter  discussed  in  the  previous  section. 
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Figure  3:  Vibration  Detector 


When  Sn{u,v )  and  S/(u,v)  are  not  known  (a  problem  often  encountered  in  pratice)  it  is 
useful  to  approximate  Eq.  (26)  by  the  relation 


F  ~  I 


l 


\H(u.v)\2 


H(u,v)  |/f(tx,u;|2  +  K 


-)G{u.v) 


(27) 


where  K  is  constant  which  is  tuned  interactively. 


5  Experiments  and  Results 

5.1  Data  Acquisition 

From  Teledvne  Brown  Engineering  Laboratory,  we  have  simulated  Aero-Optics  images  gen¬ 
erated  from  emulating  the  mixing/shear  layer  as  shown  in  Fig. 3.  These  test  images  have 
dimension  128  x  128  (6  bit  gray  levels)  and  have  two  different  set.  One  set  of  images  (128  x  128 
x  69  frames)  called  ”tare:’  were  obtained  without  mixing  CO->  gas  flow.  One  of  these  images 

is  shown  in  FigA.  The  other  set  of  images  called  ’’C02test:!  images  (128  x  12S  x  147  frames) 

Ccro-^A'vcfcr"’ 

was  obtained  with  mixing  simulated  environmental  and  the  aelUnt  COi  gas  flow.  When  the 
gas  flow  is  flowing  on  the  emulating  mixing/shear  layer,  the  impact  of  flow  may  generates  the 
vibration  and  causes  the  motion  blur.  Hence  we  detour  the  laser  beam  from  the  window  to 


Figure  4:  Pseudo  Colored  Tare  Image  (Unblurred  Image)  at  First  Frame 

Bitwise  Mask 


0 

n 

□ 

~n 

Ll 

TfftTffT 

Logical  AND  Operation 
Figure  5:  Bitwise  Logical  AND  operator 

avoid  vibration  effects.  Since  these  test  images  contain  6  bit  gray  levels,  we  took  logical  AND 
operation  with  the  bitwise  mask  as  shown  in  Fig.b.  One  of  these  "Co2test"  images  is  shown 
in  Fig.b.  The  image  shown  in  Fig. 6.  has  a  very  low  contrast. 

Image  contrast  is  related  to  the  range  of  gray  levels  in  an  image;  the  greater  the  range,  the 
greater  the  contrast  and  vice  versa.  Contrast.  C,  may  be  defined  numerically  in  the  following 


wav. 


C  =  GLmax  -  GL 


(28) 


where  GLmax  and  GLmtA  are  the  maximum  and  minimum  gray  levels  in  the  image.  Normally 


Figure  6:  C02test  Image  (Blurred  Image  by  C02  gas  at  First  Frame 

GLmax  =  -55  and  GLmtn  =  0.  Thus  image  contrast  C  is  255  (S  bit  gray  levels)  for  the  most  of 
images.  In  fact.  :,Co2test"  images  have  a  GLmax  =  10  and  GLmtn  =  0.  Since  human  eyes  can 
not  have  a  visual  discrimination  for  this  low  contrast  image  as  shown  in  Fig. 6,  one  can  not  see 
any  object  in  this  image.  Hence  contrast  manipulation  is  needed  to  increse  the  contrast  of  a 
displayed  ”C02test:’  image  by  expanding  the  original  gray  level  range  to  fill  the  dynamic  range 
of  the  display  device. 

5.2  Contrast  Manipulation 

Contrast  manipulation  is  a  pixel-bv-pixel  radiometric  transformation  that  is  designed  to  en¬ 
hance  visual  discrimination  of  low  contrast  image  features.  Each  pixel’s  gray  level  is  changed 
by  the  specified  transformation,  without  regard  for  neighboring  pixel  gray  levels.  The  contrast 
modification  performed  in  this  expriment  is  linear  radiometric  transformations.  This  simple 
linear  transformation  is  routinely  used  to  increase  the  contrast  of  a  displayed  image  by  expand- 
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New 

GLmax 


Here  H(u,v)  ,  G'(u,  u),  F{u,v)  and  N(u,v )  represents  samples  of  the  frequency  response  of 
h(i.j ),  g(i.j).  and  n(i.j)  respectively.  By  ignoring  noise  term. 


H(u.v) 


G(u,v) 
F(u,v ) 


0  <  u,  v  <  127 


(30) 


5.3.2  Image  Restoration 


Since  Nfu.vj  is  12S  x  128  dimension,  we  can  not  estimate  N(u.v)  globally.  Hence  we  use 
the  approximated  Wiener  filter  by  tuning  K  value  as  shown  in  Eq.(27).  Sometimes  one 
can  obtain  estimates  of  the  noise  variance  and  possible  noise  power  spectrum  Sn{u,v)  from 
relatively  smooth  local  regions  in  the  blurred  image  g{ij).  By  selecting  32  x  32  smoothed  local 
region,  one  can  estimate  noisy  variance  cr  by  the  following  equation. 

<7=  =  -  £{»(  U)))2  .  0  <  /,  k  <  31  (31) 

*’n  1  (=0  k—Q 

As  suggested  by  i4j.  one  can  select  the  tuning  factor  K  =  2<72.  After  Wiener  filtered  according 
to  Eq.  (27)  .  inverse  discrete  fourier  transformed,  and  linear  contrasted  restored  image  f(x.y) 
is  displayed  Fig. 9. 


For  the  mulchannel- image  case  of  interest,  let  F(u,  u),  G,-(u.  v),  and  H,(u,  v)  be  the  Fourier 
transformed  images  of  the  unblurred,  the  ith  blurred  image,  and  the  fth  point  spread  function, 
respectively.  According  to  [3],  the  Multichannel  Wiener  filtered  image  estimate  is 

(32) 


^)  =  7J»^(  ...i 


(Ei,W(“,WI3  +  A'.„s) 
where  i  =  1,2 .  N  Channels  and  Kavt<i  =  ~  with  crt-  estimated  noise  variance  of  each  i 


'aveg  —  y 

channel  imase  by  selecting  32  x  32  local  window. 
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Having  more  than  one  blurred  image  of  a  common  object  provides  information  that  can 
be  used  to  advantage  in  the  restoration  process.  One  or  the  advantages  is  in  the  area  of 
noise  reduction  through  statistical  independence  between  image  noise  fields.  The  second  is 
the  potential  elimination  of  zeros  in  the  denominator  of  the  restoration  filter  [Eq.  (27)]  by 
proper  selection  of  the  Hi(u,v)' s  or  by  requiring  just  one  Hi(u.v)  to  have  noncommon  zeros 
with  ail  other  Hi{u,v)' s.  Fig.  10  shows  the  Multichannel  Wiener  filter  restored  and  contrast 
wretched  images. 

5.4  Image  Restoration  by  Local  Image  Stastistics 

5.4.1  Blur  Identification 


This  blur  identification  technique  is  developed  by  [2],  The  degraded  image  g(i.j)  is  divided 
into  subimage  )  which  may  overlap.  If  one  assumes  that  the  extent  of  the  PSF  h(i,j) 
is  small  compared  to  the  extent  of  the  subimage  gi(i.j)  and  ignoring  the  noise  term.  then, 
approximately. 

G',-(u.  v)  H(u.v)F,(u.  v)  (33) 


According  to  (2],  the  magnitude  of  H[u,v)  is 


,ff,.  rlP  _  iES'iau.** 


(34) 


where  i  is  each  subimage  index  and  .V  is  the  total  number  of  subimages.  In  our  experiments, 
we  select  32  x  32  subimage  with  0  and  50  %  overlapping  neighborhood  subimage  respectively. 
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f 


Figure  11:  Estimated  a)  magnitude  and  b)  phase  spectrum 
The  phase  of  H(u,v)  is 

9h(u  +  Su.v  +  Sv )  =  9h(u.  y) 

+  (^F(ui  *>)  -  0Fi(u  +  Su.v  +  6u)]awe5 

Ff5r.ll  shows  the  estimated  32  x  32  magnitude  and  phase  spectrum  of  H(u.v). 
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5.4.2  Image  Restoration 


After  identifying  the  magnitude  and  phase  spectrum  of  H(u,v),  we  can  apply  32  x  32  dimension 
Wiener  filter  H'32x32(u5u) 


^32  X32(U>U)  — 


Hm(u,v) 


(35) 


lff(u,u)|2  -f  '± 

^n32x32  (^»v) 

where  Sj22xi2{u.  v)  denote  power  spectral  densities  of  the  undegraded  image  and  (u,u) 

denote  power  spectral  densities  of  noise  by  selecting  local  smoothed  32  x  32  pixels  region. 
Estimating  power  spectral  densities  by  using  Fast  Fourier  Transform  is  well  documented  in  [9]. 
In  our  experiments,  these  power  spectral  densities  are  estimated  by  the  following  Periogram 
method. 

1  K  |  „\\2 

(36) 


°/33x33  ~  ,V  Z- 


Nfr[  32x32 

where  N  is  the  total  number  of  subimages  and  F;(u,  u)  is  Fourier  transformed  tth  subimage.  Af¬ 


ter  selecting  smoothed  32  x  32  pixels  region  in  the  blurred  image, 
by 

liY(u.u)l2 
n5iM<3J'  32x32 

where  N(u.v)  is  Fourier  transformed  noise  region. 

The  estimated  image  Ft(u,u)  is  given  by 


'"33X33 


(u.v)  can  be  estimated 


(37) 


Fe(u,u)  =  W32x32(u,v)G{u,v) 


(38) 


where  G(u,v)  is  Fourier  transformed  blurred  image.  But  the  dimension  of  ir32X32(u. v)  >s  32  x 
32  and  that  of  G(u.v)  is  12S  x  12S.  Hence  we  need  to  interpolate  (or  zoom)  H?2x32(“*t;)  into 


128  x  128  pixels  dimension.  Several  image  interpolation  techniques  are  discussed  in  [6].  In  our 
experiments.  bilinear  interpolation  technique  is  used.  This  algorithm  uses  the  four  input  pixels 
surrounding  the  point  (t.jj  to  estimate  the  output  pixel.  After  interpolate  W„x:a(u.v),  one 

can  take  Inverse  Fourier  transform  to  get  restored  image.  The  restored  and  contrast  stretched 
images  are  shown  in  Fig.  12  and  Fig.  13. 


Figure  13.  (c) :  Restored  Image  (Pseudocolored)  by  50  Percent 

overlapping  subimage  at  69th  Frame. 


5.5  Conclusion 


In  this  report,  we  can  identity  the  mixing/shear  layer  blurring  process  nonparametrically  by 
estimating  the  magnitude  and  phase  spectrum  of  H{u,  v).  By  applying  Wiener  filter  using 
global  image  statistics,  we  have  restored  image  from  the  seveley  degraded  image  (See  Fig. 6). 
In  this  restored  image  ,  the  original  target  is  clearly  visible  on  the  center.  But  it  has  a  false  slit 
pattern  along  the  vertical  axis  and  has  one  CO 2  gas  pattern.  This  false  slit  pattern  might  be 
i'ome  from  utilizing  global  image  statistics  when  one  estimate  Wiener  filter.  In  the  restored 
image  as  shown  in  Fig. 12  and  Fig.  13.  several  distinct  CO2  gas  patterns  or  objects  are  visible 
and  bright  spots  related  with  original  target  are  discerned  around  the  center. 
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CONTENTS 


a)  Task  Objective:  The  overall  purpose  of  the  research  is  to  create  an  algorithm  which 
can  generate  optical  or  infrared  images  from  a  given  source  when  the  image  passes  through 
a  thin  layer  which  is  highly  turbulent.  The  turbulent  layer  is  assumed  to  be  induced  by 
coolant  passing  over  a  hypervelocity  optical  window.  The  algorithm  should  be  capable  of 
studying  image  effects  when  the  turbulence  is  inhomogeneous  and  anisotropic.  The 
algorithm  should  also  be  capable  of  studying  image  effects  as  arise  when  the  image  passes 
through  a  bow  shock  region  which  precedes  the  turbulence  layer  and  also  study  image 
effects  as  they  occur  when  the  image  passes  through  a  heated  aero-optical  window  which 
occurs  after  the  image  has  passed  through  the  turbulence.  The  algorithm  wiil  be  very 
useful  for  gauging  how  badly  different  hyperveiodty  aero-optical  scenarios  wiil  effect  image 
quality  of  different  optical  sensors. 

b)  Technical  Problems:  The  level  of  effort  is  mainly  involved  with  theoretical  optical 
analysis  and  numerical  analysis.  The  main  technical  problems  have  concerned  lack  of 
adequate  computer  software  on  the  UAH  ECE  ARDENT  computer  to  perform  analysis. 

c)  General  Methodology:  The  general  methodology  has  been  a  theoretical  and 
numerical  optical  analysis  of  image  propagation  in  an  inhomogeneous  turbulence  layer. 
The  details  of  the  analysis  are  given  in  the  report. 

d)  Technical  Results:  See  attached  report. 

e)  Important  Findings:  See  the  conclusion  of  the  attached  report. 


abstract 


The  propagation  of  the  Mutual  Coherence  Function  is  governed  by  a  Partial 
Differential  Equation.  This  equation  can  be  solved  by  standard  Finite  Difference 
Techniques.  For  the  case  of  propagation  through  homogeneous  turbulence,  it  is  also 
possible  to  obtain  a  closed  form  integrated  solution.  We  show  that  it  is  necessary  to  split 
the  Mutual  Coherence  Function  into  a  delta  and  non-^elta  component.  We  obtain  a 

differential  equation  for  the  non^elta  component  and  compare  the  results  obtained  by  an 
iterative  solution  and  integrated  solution. 


NOTATION 


E(Z,  Xj) 

rs(z,  Xr  Xj) 


rs(z.  Xc  xd> 


V2'  X,.,  k) 


♦„(*) 

»(xd> 


=  V=r 

=  electric  field  at  a  point  (Z,  X^) 

=  <E(Z,  Xx)  E*(Z,  X2)> 

=  ensemble  average  if  the  product  of  the  electric  fields  at  E(Z,  X^) 
and  E(Z,  X2) 

=  Spatial  Domain  mutual  coherence  function  at  (Z,  Xc,  X^) 

Xt  +  X, 

“ - 2 — 1 

X1  "  X2 

= - 3 - 

=  Fourier  Domain  MCF 

CD 

=  s/rs(z'xcxd)“PHvxd)dXd 

=  wavelength  of  radiation  (m-1) 

=  2ir/\ 

=  inner  scale  of  turbulence  (m-1) 

=  outer  scale  of  turbulence  (m~*) 

=  structure  constant 
= 

=  spatial  spectrum  of  the  turbulence 

CD 

=  /  <J>n(«)  exp(-j«Xd)  d* 
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1.  INTRODUCTION 

The  image  of  an  object,  suffers  an  inevitable  degradation,  when  it  propagates 
through  a  turbulent  medium.  Such  a  medium  is  characterized  by  a  randomly  varying 
index  of  refraction.  In  the  spectral  domain,  this  varying  index  of  refraction  may  be 
modeled  by  the  Von  Karman  spectrum,  $n(«). 

The  Mutual  Coherence  Function  (MCF),  is  an  important  tool  with  which  to  analyze 
the  distortion  suffered  during  the  process  of  image  formation.  The  propagation  of  the  MCF 
is  governed  by  a  first  order  hyperbolic  partial  differential  equation,  which  is  applicable 
universally  to  the  study  of  free  space  propagation  and  propagation  through  turbulence. 
The  differential  equation  describes  how  the  MCF  is  dependent  on  the  spectrum  of  the 
refractive  index. 

In  this  report  we  are  interested  in  studying  the  problem  of  what  happens  to  an 
image  when  it  passes  through  a  thin  layer  of  material  which  is  highly  turbulent.  For 
present  application  the  thin  layer  of  material  is  formed  by  coolant  passing  over  an 
aero-optical  window.  In  our  formulation,  a  coherent  point  source  is  situated  at  a  distance 
L_,  from  a  layer  of  turbulence  of  width  L  .  We  analyze  the  image  formed  at  the  piane  Lg  + 
L,  by  first  computing  the  MCF  in  this  piane.  Once  the  image  for  a  point  source  has  been 
obtained,  the  image  for  any  general  object  intensity  distribution  can  be  obtained  by 
superposition  from  the  point  source  solution. 

At  the  present  time,  a  closed  form  integrated  solution  is  already  available  to 
compute  the  MCF  for  simple  geometries  like  this.  In  this  report  we  develop  an  alternative 
iterative  solution  for  the  MCF  PDE,  and  compare  the  integrated  and  iterative  solutions. 
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Fig.  1  Propagation  through  Homogeneous  Turbulence. 
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Propagation  through  Inhomogeneous  Turbulence. 


The  motivation  for  developing  an  iterative  solution  is  so  that  in  future  analysis, 
layers  of  turbulence  which  are  transversely  inhomogeneous  can  also  be  studied  with  the  aid 
of  the  iterative  MCF  formulation.  The  presently  available  integrated  solution  cannot  be 
used  for  this  purpose. 

We  show  that  the  MCF  can  be  expressed  as  the  sum  of  a  delta  component  and  a 
non-delta  component.  This  decomposition  is  shown  to  be  an  essential  step  for  the 
practical  numerical  implementation  of  the  iterative  and  the  integrated  solution  methods. 
For  a  coherent  point  source,  only  the  delta  component  is  present  for  free  space  propagation. 
As  the  wave  travels  through  turbulence,  the  non-delta  component  gradually  gains 
strength,  while  the  delta  component  wanes. 

Further,  it  is  shewn  in  the  report  that  the  deita  component  represents  the  true 
image  (without  turbulence)  of  the  point  source,  while  the  non-delta  component  contributes 
to  the  blur.  Thus,  as  the  wave  propagates  through  the  layer,  the  strength  of  the  image 
point  decreases,  while  the  blur  increases,  due  to  the  gain  in  strength  of  the  non-delta 
component.  Suggestions  of  how  to  apply  the  present  theory  to  transversely  inhomogeneous 
layers  is  given. 

2.  MATHEMATICAL  PRELIMINARIES 

A.  The  Differential  Equation  for  the  Propagation  of  the  MCF: 

As  shown  in  Fig.  1,  a  point  source  of  light  is  situated  at  a  distance  L  from  a  layer 
of  turbulence.  The  direction  of  propagation  of  the  light  is  orthogonal  to  the  turbulent 
layer.  Sensors  are  located  at  a  point  Z  =  L  in  the  turbulent  layer.  We  wish  to  find  the 
MCF  in  the  plane  of  the  sensors.  The  turbulence  is  governed  by  the  Von  Karman 
spectrum,  <L(k)  (1,  2j.  The  MCF  in  the  spatial  domain  ?  ,  is  governed  by  the  following 
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partial  differential  equation  (1,  2]. 
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2jklz+  (fflcfaxlJ  -irWo)~a(xi  -x2>]}  rs(z,  xr  Xj)  =  o 


Where,  f S(Z,  Xj,  Xj)  _  <E(Z,  Xj)  E  (Z,  Xj)>  and  other  variables  are  defined  in  the 
-Notation  table.  Now  transforming  to  the  (Xc,  Xd)  coordinate  system,  where, 


Xj  +  X,  X,  -  X, 

Ac  2  an°  xd  = - 2 - we  find 


(  2  3 

l2jk  th  +  5  +  jfw°)  -  *(xd)]}  q(z,  xc,  xd)  =  0 .  (2) 

To  proceed  further  we  wiU  express  Eq.  (2)  in  the  spectral  domain  through  the  use  of 
the  following  spatial  Fourier  transform: 


rs(z>  xc>  xd)  =  /  r^Z,  Xc,  Kc)  exp(-j«cXd]  d/c. 


a(xd)  =  /  $n{K)  exp(-j«X.)  d* 


2/2 


♦aw = °-033  cJ(«2 + Vu/6  e~*  fK- 

Lo 

m  T 


After  substituting  (3)  in  (2),  we  apply  the  well  known  result,  that  the  product  of  spatial 
domain  terms,  TJZ,  Xc,  Xd)  and  a(Xd)  is  equivalent  to  the  convolution  of  the  frequency 
domain  terms,  Tj(Z,  Xc,  «c)  and  *  (<t)  and  obtain: 


(5a) 


This  equation  describes  the  propagation  of  the  MCF  through  homogeneous  turbulence. 

Since  the  focus  of  the  report  is  on  the  study  of  propagation  through  inhomogeneous 
turbulence  as  weil  as  homogeneous  turbulence,  we  will  now  modify  Eq.  (5a)  to  account  for 
the  more  general  inhomogeneous  case.  Fig.  2  shows  a  more  general  case,  where  the 
turbulence  layer  is  wedge  shaped.  To  solve  (5a),  for  turbulence  of  inhomogeneous 
geometry,  we  multiply  the  last  term  of  Eq.  (5a)  with  a  function  f(Z,  Xc): 


dr  3 

_2jk  ^ -  2j«c  -fij- -  jy-  f(Z,  Xc)  |a(0)rf -  f  I^Z,  Xc,  «^)<|)n(«c  -  K'c)  d/c^.  j  =  0  (5b) 


We  have  used  the  Von  Karman  spectrum  to  describe  <J>  (/c)  and  we  model  f(Z,  X  )  by  a 

la  U 

polynomial  function,  which  takes  a  value  zero  in  free  space  and  gradually  attains  a  value  of 
1,  as  it  enters  deeper  into  the  turbulence. 


3.  The  Integrated  Solution  for  the  Homogeneous  Case  [1]  [2]: 

An  integrated  solution  r^Z,  XQ)  *c)  of  the  (5),  with  f(Z,  Xc)  =  L  (homogeneous 
case)  has  been  shown  to  be  (1,  2]: 
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Although  mathematically  correct,  computationally,  (6)  is  not  easy  to  handle,  because  it 
contains  a  delta  component,  and  hence  the  integral  in  (6)  does  not  converge  in  the  ordinary 
sense.  To  see  this,  we  note  that  in  general 


Therefore,  we  must  find  a  way  to  separate  the  delta  and  non-delta  components  in 
fj(Z,  X£,  /cc)  Eas.  (6)  and  (7)  may  aiso  be  written  in  the  form: 


K  2 

r j(Z,  Xc>  kc)  =  ^  |  f  exp(j*d  ^  Z  -  X£  ]  exp  [-  ~~  a(0)(Z  -  Ls)] 


— B 


[exp[£/a(: 


ka  /I 


dZ 


')  -  iK( 


*  /t  2 

ifife  |  /exp[jxd  y  Z  -  Xc] ]  exp  [- 1 -  a(0)(Z  -  Lg)] J 


(8) 


The  second  term  in  Eq.  (8)  is  a  weil  known  representation  of  a  delta  function.  Using  this 
information  we  may  write  as  a  sum  of  a  non-delta  component  Tj  and  a  delta  component 
T. 


iyz,  xc,  *c)  =  r'(z,  xc,  «c)  +  t(z,  xc,  kc) 


(9) 


where 
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Fig.  3 


Plot  of  a(Xd)  versus  Xd< 


(10) 


rJ(Z,  alig¬ 


ns  {  /expu«d  [j£  Z  -  Xc]  1  exp  t- 1-  a(0)(Z  -  Ls)] 


T(Z,  Xc,  *c)  =  £  expj-  a(0)(Z  -  L,)] 


z-x. 


(11) 


As  opposed  to  Eq.  (5),  the  Fourier  Transform  in  Eq.  (10)  does  converge.  Consider  the  last 
term  of  the  product. 


2  ?  m 


K  (exP  T~  J  a  I  —  ^1  ~  exP  (®)  “ l  —  0 


Since  a(X)  approaches  zero  as  X  approaches  infinity  (See  Fig.  3),  the  integral  in  Eq.  (10) 
does  converge. 

A  further  property  of  the  integrated  solution  is  that  we  need  to  carry  out  the 
integral  in  Eq.  (10)  only  for  XQ  =  0.  This  results  because  of  the  form  of  k  and  Xc  in  the 
argument  of  the  first  exponential  of  Eq.  (10).  If  r^Q  is  solution  for  r  £  at  X  =  0.  then 
when  X  t  0  is  given  by 


rr0<z'  *c>  =  ri(z'  °- 

rf(Z,  Xc.  *c)  =  r'Io(Z,  Kc  - !  Xc).  (12) 

Next  we  develop  an  iterative  solution  for  T^.  It  is  one  of  the  goals  of  our  work  to  compare 
the  iterative  solution,  with  the  integrated  solution,  given  by  (10). 
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C.  The  Iterative  Solution 


To  find  the  iterative  solution  for  the  homogeneous  case,  we  set  f(Z,  X  )  =  1  and 

apply  the  PDE  given  by  (5),  to  each  of  the  components  of  I^Z,  Xc,  kq).  Substituting 
=  r  f  t  T  into  Eq.  5a  we  have: 


jZ(rt  +  T)+j£|3r(r' +t) 


+ 


2  ® 

f{a(0)(r;+T)-/[r;(Z,Xc,^)+T(Z,Xc,K')]^(VK.)d«; 


=  0 


(13) 


Rearranging  terms. 


( 5T  *c  5T  T  * 

W + r  <irc  +  r  H°)  t  -/t(z,xc,*'  )<U' 


arf  1 2  r 

3Z-+r33q+r[1(°)rr/rf(z.xc.''c»n(V'‘c)d'tc]} =  °  (u) 


Denoting  the  terms  in  the  first  set  of  braces  by  Q(Z,  X£,  *c)  and  rearranging  (14): 


dr  ^  ^2/ 

32T+ r  r{sWrr/rf(z.xc''!i)+n(V'!Pd^}  =  -Q(Z.  xc,  *( 


)  (15) 


and 


Q(Z'Vc>  -  3T +  r  JTC  +  r [a(°)T~/T(z, V'  »„( V"c)d^ 


(16) 
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We  use  Eq.  (11)  for  T  and  obtain  and 

c 


-  -  f  a(0)  exp  [-  £  a(0)(Z  -  Ls)]  S  (£  Z  -  X 

Kr  T  1  f*r 

+  -^exp  -^-a(0)(Z - Lg)|  S'  |p- 


k  k  r  Ko  i  i 

r  =  _  ri exp  ["  aMz  ~  Ls)J 6>  [r  z  "  xc 


Z"XcJ 

(17) 

(rz_xc 

(18) 

Using  (17)  and  (18)  in  (16),  we  2nd  that  the  delta  functions  drop  out  and  that  0  is  given 
by 


Q  = 


(19) 


Again  using  the  expression 


JU  *  ••  Ay 


2  2 

Q  (2.  Xc,  *c)  -  f  exp  {-  ]j-  al  0)(Z  -  y}  -  \  Xc) 


(20) 


The  nonhomogeneous  right  hand  side  of  Eq.  (15)  which  results  from  integration  over 
the  T  delta  function  can  be  thought  of  as  a  source  which  excites  as  the  wave  propagates 
in  the  medium.  Because  T  of  Eq.  (11)  decreases  with  increasing  Z,  one  can  think  of  energy 

as  being  transfered  from  the  T  delta  function  (coherent  energy)  to  the  incoherent  part  of 

» 

the  wave  T  ^  (incoherent  part  of  the  wave).  It  is  interesting  that  the  spectrum  of  the 

/ 

turbulence  itself  determines  how  the  incoherent  energy  T  ^  is  generated.  Now  substituting 
(20)  into  (15)  gives 
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(21) 


a> 


dr  [  KcdF[  ,2  , 

tt+ r °rc  +  r |a<0>rf -  /rf(z. xc,  *-)  v-c-^)  a*d 

-« 

2  2 

=  iz  “p  {-  r  a(°)(z  -  Ls)}  ♦>«.  -  \  xc) 

We  now  discrete  this  equation  at  the  point  (Zn+1>  Xem,  v)  and  we  define 


AZ  =  Z_  ,  t  - Z_  AX  =  X  ,  _  x  Alt  -  K 

r*  i  1  -»■ -  ^  4  ,  —  ft 


n+1  a 


c  '  cm+l  *  cm  a*c  “  *c/+i  ~  *c^ 


(22) 


The  discretised  system  now  becomes: 


The  initial  condition  is: 


rf(V\.*c>  =  °  (24) 

The  system  (23)  is  a  backward  difference  scheme  for  approximating  PDE’s.  The 
convolution,  when  discretised  leads  to  the  well-known  Teplitz  matrix,  whose  elements  are 

various  values  of  the  function  a(X4).  This  matrix  remains  the  same  throughout  the 
computation  and  needs  to  be  inverted  only  once. 
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3.  EXAMPLES  AND  RESULTS 

In  this  section,  we  present  a  few  examples  and  resuits.  We  solve  a  homogeneous 
problem,  typified  by  Fig.  2,  and  compare  the  results  obtained  from  the  iterative  solution, 
which  is  set  up  in  equation  (23),  with  the  results  obtained  from  the  integrated  solution  in 
Eq.  (10). 

For  all  the  examples  we  use  the  following  values  for  the  parameters: 

Ls  =  1000  m 

A  =5.5*  10"7  m  (5500*  A) 

C  =  0.2  m-1/3 
n 

£  =  10“5  m 

o 

L  =  0.001  m 
o 

AX.  =  10-3  m 
c 

AZ  =  10-4  m 
Ak,  =  50  m-1 

Fig.  4  shows  how  the  strength  of  the  delta  function  T(£,  0,  0)  drops  as  it  propagates 
through  the  turbulent  layer.  Z  is  the  distance  in  the  turbulent  layer  and  Z  =  L  fZ. 

Fig.  5  is  a  plot  of  the  non-delta  component,  Tf(2,  0,  0)  (Z  =  L  t  Z,  X  =  0.  k  = 

I  5  L  V* 

/ 

0).  There  are  two  completing  processes:  continually  gains  energy  from  the  delta 

component  T(Z,  Xc,  *c),  but  also  loses  energy  as  it  propagates  through  the  turbulence. 
Since,  we  have  seen  in  Fig.  4,  that  the  strength  of  T(Z,  X  ,  k  )  continually  decreases,  there 

comes  a  point  where  the  loss  of  energy  to  the  turbulence  cannot  be  offset  by  a  gain  of 

/ 

energy  from  T(Z,  XQ,  kq)  and  Tj  starts  to  decrease. 

Fig.  6  is  a  3-D  plot  of  0,  *c)  from  £  =  0toZ  =  3*  10~3  mm. 

Fig.  7  is  a  3-D  plot  of  Tj(2,  10,  *c)  and  with  AXC  =  0.001  metres,  this  shows  how 
changes  at  off  center  point,  meters  away  from  the  origin. 
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Comparison  of  the  integrated  and  iterative  solutions  at  a  distance 
i  *  3  *  10  meters  in  the  layer.  The  continuous  curve  is  the  iterative 
solution,  while  the  asterisks  represent  the  integrated  solution. 


o 

Fig.  8  is  a  comparison  of  the  iterative  and  integrated  solution  for  Z  =  3  x  iq  J  m, 

=  0.  The  continuous  curve  is  the  iterative  solution,  while  the  asterisks  represent  the 
integrated  solution; 

4.  CONCLUSIONS 

We  have  solved  the  differential  equation  governing  the  propagation  of  the  Mutual 
Coherence  Function  in  the  frequency  domain,  by  two  techniques:  an  iterative  method  and 
a  closed  form  integrated  method.  The  results  from  both  these  methods  agree  extremely 
well. 

It  was  necessary  to  split  the  MCF  into  two  components:  a  delta  component  and  a 
non-delta  component.  In  free  space,  only  the  delta  component  is  present  and  a  point 

source  images  to  a  point,  if  diffraction  at  a  lens  aperture  is  ignored. 

/ 

In  a  turbulent  medium,  the  non-delta  component  is  a  significant  factor.  It 

grows  in  strength  by  absorbing  energy  from  the  delta  component,  until  it  reaches  a 

maximum.  Simultaneously,  it  also  loses  energy  to  the  medium  by  scattering.  After 

/ 

reaching  a  maximum  gradually  falls  off. 

The  delta  component  will  stiil  image  to  a  point,  with  greatly  diminished  strength, 

.  / 

while  the  Tf  will  contribute  to  fuzziness  in  the  image. 

5.  FUTURE  WORK 

We  intend  to  study  the  propagation  of  the  MCF,  for  inhomogeneous  turbulence. 
We  will  then  construct  an  image  of  a  point  source,  after  propagation  of  the  wave  through 
the  inhomogeneous  turbulence. 

We  have  concentrated  our  efforts,  so  far  on  two  dimensional  models.  The  next  step 
will  be  to  analyze  the  more  difficult  three  dimensional  case. 
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Fig.  1 

Fig.  2. 

Fig.  3 

Fig.  4 

Fig.  5 

Fig.  6 

Fig.  7 

Fig.  8 


Figure  Captions 

Propagation  in  Homogeneous  Turbulence. 

Propagation  in  Inhomogeneous  Turbulence. 

Plot  of  a(X(j)  versus  X^. 

T(Z)  at  X  =  0,  k  =  0  versus  Z,  where  Z  =  distance  in  the  turbulent  layer 
and  Z  =  t  Z. 

Ff(Z)  at  X£  =  0,  /c  =  0  versus  Z.  Z  =  Lg  +  Z. 

A  3-D  plot  of  rf(Z,  «c)  at  Xfi  =  0  meters. 

A  3-D  plot  ofrf(Z,  «c)  at  X£  =  0.01  meters. 

Comparison  of  the  integrated  and  iterative  solutions  at  a  distance 
—3 

Z  =  3  *  10  meters  in  the  layer.  The  continuous  curve  is  the  iterative 
solution,  while  the  asterisks  represent  the  integrated  solution. 
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